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A new algorithm for the symboHc computation of polynomial conserved densities for 
systems of nonlinear evolution equations is presented. The algorithm is implemented 
^ • in Mathematica. The program condens.m automatically carries out the lengthy sym- 

OP ' bolic computations for the construction of conserved densities. The code is tested on 

s»J several well-known partial differential equations from soliton theory. For systems with 

>«.i^ , parameters, condens.m can be used to determine the conditions on these parameters 

C_N ' SO that a sequence of conserved densities might exist. The existence of a large number 

^a_^ ' of conservation laws is a predictor for integrability of the system. 

!>■ 



1. Introduction 

^ , An indication that certain evolution equations might have remarkable mathematical 

properties came with the discovery of an infinite number of conservation laws for the 

- T— , , Korteweg-de Vries (KdV) equation, ut+uux+u^x = 0. The conserved quantities u and u^, 

^\^ ' corresponding to conservation of momentum and energy, respectively, were long known, 

3 , and Whitham (1974) had found a third one, m'^ — 3u^, which corresponds to Boussinesq's 

moment of instability. Zabusky and Kruskal found a fourth and fifth. However, the search 

for additional conserved densities was halted due to a mistake in their computations. 

Miura eventually continued the search and, beyond the missing sixth, found an additional 

three conserved densities (Newell, 1983). It became clear that the KdV equation had an 

infinite sequence of conservation laws, later proven to be true. 

The existence of an infinity of conserved densities was an important link in the dis- 
covery of other special properties of the KdV equation (Zakharov, 1990). It led, for 
example, to the construction of the Miura transformation, which connects the solutions 
of the KdV and modified KdV equations. Consequently, the famous Lax pair was found, 
which associates a couple of linear equations to the KdV equation. From that, the inverse 
scattering technique (1ST) for direct linearization of integrable equations was developed, 
and it was then shown that the KdV equation, and many other integrable equations, 
admit bi-Hamiltonian structures. 
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There are several motives to construct conserved densities of partial differential equa- 
tions (PDEs) explicitly. The first few conservation laws have a physical interpretation. 
Additional ones may facilitate the study of both quantitative and qualitative properties 
of solutions (Scott et ai, 1973). Furthermore, the existence of a sequence of conserved 
densities (perhaps with gaps) predicts integrability. Yet, the non-existence of conserved 
quantities does not preclude integrability. There are indeed equations, viz. dissipative 
equations, with only one conserved density, which can be directly integrated. The most 
notable is the Burgers equation, which can be transformed into the linear heat equation 
via the Cole-Hopf transformation (Zakharov, 1990). 

Yet another compelling argument to explicitly construct conserved densities relates to 
the numerical solution of PDEs. It is desirable that the semi-discretization conserves the 
discrete analogues of the continuous conserved quantities. In particular, the conserva- 
tion of a positive definite quadratic quantity may prevent the occurrence of nonlinear 
instabilities in the numerical scheme. The use of conservation laws in solving the Boussi- 
nesq equation numerically has been illustrated in Hickernell (1983). Sanz-Serna (1982) 
describes a scheme for the integration in time of PDEs, which is explicit and capable of 
conserving discretized quadratic functionals. Since u and u^ are conserved densities for 
the KdV equation, a discrete scheme should have conservation of momentum, ^ t/", and 

energy ^ [[/"] . In Sanz-Serna (1982), an explicit self-adaptive conservative scheme with 
conservation of energy and momentum is given. Conservation of energy implies bound- 
edness of the solutions, and therefore obviates the occurrence of any blowup phenomena. 
For more details about numerical applications of conservation laws see LeVeque (1992). 

We present a new symbolic algorithm for computing closed-form polynomial-type con- 
servation laws for systems of nonlinear evolution equations. Our algorithm is also ap- 
plicable to wave equations, viz. the Boussinesq equation, provided that the PDE (of 
higher-order in time) can be written as a system of evolution equations (first order in 
time). In contrast to fairly complicated algorithms designed by Bocharov (1991), Gerdt 
and Zharkov (1990), and Gerdt (1993), we introduce an algorithm that is based in part 
on ideas presented in Hereman and Zhuang (1995), Ito (1994), and Ito and Kako (1985), 
Kruskal et al. (1970), Kodama (1985), Miura et al. (1968), Verheest and Hereman (1995), 
and Willox et al. (1995). Our algorithm has the advantage that it is fairly straightforward 
to implement in any symbolic language. We also present a software package condens.m, 
written in Mathematica syntax, which automates the tedious computation of closed-form 
expressions for conserved densities and fluxes. 

In Section 2, we give the definitions of conservation law, density and flux. We also 
state a theorem from calculus of variations about the Euler-Lagrangc equations, which 
will play a role in our algorithm. The remainder of the section is devoted to a detailed 
exposition of the algorithm, which was implemented in Mathematica, and successfully 
tested on many well-known evolution systems from soliton theory. 

In Section 3, we address applications of our program condens.m. In particular, we 
show how it could be used in a search for integrable fifth-order evolution equations of 
KdV type, where we retrieved all the known integrable cases. We carried out a similar 
computer search for a parameterized class of seventh-order evolution equations. More 
examples and test results are presented in Section 4. In Section 5, we describe the usage 
of our code, and indicate its capabilities and limitations. In Section 6, a comparison with 
other programs is given. Also, several ongoing projects arc briefiy addressed. 
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2. Computation of Conserved Densities 

2.1. Definitions 

For simplicity, consider a single PDE, 

A{x,t,u{x,t)) = 0, (2.1) 

where t G M denotes time, x (^ Mis the spatial variable, and u{x, t) G M is the dependent 
variable. A conservation law is of the form 

Dtp + D,J = 0, (2.2) 



which is satisfied for all solutions of (2.1). The functional p(x, t) is the conserved density, 
J{x,t) is the associated flux. Both are, in general, functions of x,t,u, and its partial 
derivatives with respect to x. Furthermore, Dt denotes the total derivative with respect 
to t; Da; the total derivative with respect to x (Ablowitz and Clarkson, 1991). Specifically, 
p is a local conserved density if p is a local functional of u, i.e. if the value of p at 
any x depends only on the values of u in an arbitrary small neighborhood of a;. If J 
is also local, then ( p.2| ) is a local conservation law. In particular, if p is a polynomial 
in u and its x derivatives, and does not depend explicitly on x or t, then p is called 



a polynomial conserved density. If J is also such a polynomial, then (2.2) is called a 



polynomial conservation law. There is a close relationship between constants of motion 



and conservation laws. Indeed, for polynomial-type p and J, integration of (2.2) yields 



P = p dx = constant, (2-3) 

provided that J vanishes at infinity. For ordinary differential equations, the P's are 
constants of motion. 

Example 2.1. The most famous evolution equation from soliton theory, the KdV equa- 
tion (Miura, 1968), 

ut + uux + u^x = 0, (2.4) 

is known to have infinitely many polynomial conservation laws. The first three are 

[u^+i^u' + u^A =0, (2.5) 

(m2)^+ Qu3_u^2^2mM2xJ =0, (2.6) 

{u^ -'iUx^)^+ {-U^ -QuUx^ + 'iu^U2x + 'iU2x'^ ~^UxUZx\ =0. (2.7) 

The first two express conservation of momentum and energy, respectively. They are easy 
to compute by hand. The third one, less obvious and requiring more work, corresponds 
to Boussinesq's moment of instability (Newell, 1983). 

Observe that the KdV equation and its densities p = u,v? and u^ — iux^ are all 
invariant under the scaling symmetry 

{x,t,u) -^ {Xx,X^t,X~^u), 
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where A is a parameter. Stated differently, u carries the weight of two derivatives with 
respect to x; denoted synibohcally by 

^^ 
dx"^' 

Scahng invariance, which is a special Lie-point symmetry, is an intrinsic property of 

many integrable nonlinear PDEs. Our algorithm exploits this idea in the construction of 

conserved densities. 

2.2. EuLER Operator 

We introduce a tool from calculus of variations: the Euler operator, which is very useful 
for testing if an expression is a total derivative (Ito and Kako, 1985; Olver, 1986), without 
having to carry out any integrations by parts. 

Theorem 2.1. /// = /(x,yi, . . . ,yj"\ . . . ,yAr, . . . ,y^^), then C^{f) = 0, if and only if 

p d („_1) (n-l)\ 

J = -^g, where g = g(x,yi,...,y{ ' , . . . ,yM, ■ ■ ■ ,yN )■ 

In this theorem, for which a proof can be found in Olver (1986, pp. 252), y — [j/i, . . . , j/at]"^, 
Cjj{f) = [Cy^ (/),..., Cy^ {f)V , = [0, . . . , 0]-^, with T for transpose, and where 
d d d (P d d"- d 

is the Euler operator (or variational derivative) . We will use this theorem in our algorithm. 

2.3. Algorithm 



We now describe our algorithm to compute polynomial conservation laws (2.2) for 
systems of nonlinear evolution equations of order n, 

Hf ^f{u{x,t),u'{x,t),...,u^"\x,t)), (2.8) 

where u = [ui, . . . , un]'^ , or, component-wise, 

u»,t+^»(uj,Uj,Uj,...,M5"^)=0, 2 ==1,2,..., TV, j = l,2,...,7V, (2.9) 

where 

''* " at ' ' ~ ■'■"" " ax" ' 

and all components of m depend on x and t. 

Our goal is to compute the densities p{u, ..., u^"^^^) and the fluxes J{u, ..., ■u'^™^)) of order 



mi and m2, respectively. There will be no restriction on the order n of the system (2. 



or on its degree of nonlinearity. But, all the components of JF must be polynomial in 
their arguments. Furthermore, we only consider systems of evolution equations in t with 
one spatial variable x. 

In our algorithm, we tacitly assume that we have an evolution equation for every 
dependent variable Ui. In cases where there are more dependent variables than equations, 
one can always add trivial evolution equations Ui^t = 0. Further on, we use the notation 
Uj^nx instead of Uj^"' because it is closer to the notation in our code. 
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Step 1 Determine the weights of variables and parameters 

We define the weight of a variable as the number of partial derivatives with respect 
to X the variable carries, and the rank of a term as the total weight in terms of partial 
derivatives with respect to x (Kruskal et ai, 1970; Miura et ai, 1968). The rank is 
assumed to be non neg ative and rational. 

For the system ( |2.9|) , we first try to determine the weights (scaling properties) of all 
the variables. We assume that all terms in a particular equatio n ha ve the same rank. 
We call this property uniformity in rank. Different equations in (p.9[ ) can have different 
ranks. Having defined the weight of a variable in terms of x— derivatives, we (can) pick 
w{-^) = 1, . . . , w{-§^) = n, where w returns the weight of its argument. 

For simple systems, in particular those with uniform rank, only the variables Ui and 
^ have weights. However, to be able to handle more general systems, we allow for 
constant parameters to be introduced and also for these parameters to carry weights. 
The trick of introducing parameters with weights allows one to handle equations with- 
out uniform rank. Let us assume that there are P such parameters in the system, de- 
noted by Pi, i = 1,2, ...,P. Thus, the extended list of variables that carry weights is 

{^,Ml,M2, • • ■ , UN,Pl,P2, ■ ■ ■ ,pp}- 

Weights must be nonnegative and rational, except for w(^), which may be positive, 
zero, or negative. More precisely, the weight of at least one Ui must be positive; the 
weights of the parameters pi must be nonnegative. We proceed with the following steps: 



(a) Take the i equation in (2.9). Denote the number of terms in the equation by Ki. 

(b) Compute the rank ri^k of the fc*'' term in the i*'' equation as follows: 

rj,fc = d{x) +d{t) w{ — ) + ^g{uj) w{uj) +^g{Pj) w{pj), k ^ 1,2, . . . ,K^, 

where g returns the degree of nonlinearity of its argument, d returns the number 
of partial derivatives with respect to its argument. For evolution equations d{t) is 
either zero or one. 

(c) Assuming uniformity in rank in the z"* equation, form the linear system 

A = {ris = ri^2 = • • ■ = n^K,}- 



N 



(d) Repeat steps (a) through (c) for all of the equations in (2.9). 

(e) Gather the equations Ai to form the global linear system A — M ^i- 

(f) Solve A for the N + P + 1 unknowns w{uj), w{pj) and w{-§^). 
If the solution of the system A still has free, consider two cases: 

1) If two or more weights are undetermined, prompt the user to enter choices. 

2) If only one weight is free, say w{uk), take the equations obtained in (f), set their 
left hand sides equal to one, and solve them piece by piece for w(uk)- Include the 
choice w{uk) = 1. For all choices for w(uk) test: (i) if w{uk) is negative, increment 
it until it is positive; (ii) reject w{uk) if any other weight is negative. Continue 
with the smallest integer value for w{uk) if present, else take the smallest fractional 
value. This produces at most one positive value for the free weight, out of possibly 
infinitely many choices. If the algorithm fails, prompt the user to enter a value for 
w{uk). 
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Step 2 Construct the form of the density 

The second step involves the construction of the polynomial density with a prescribed 
rank R. All the terms in the density p must have that same rank R. Since we may 
introduce parameters with weights, the fact that the density will be a sum of monomials 
of uniform rank does not necessarily imply that the density must be uniform in rank 
with respect to the depen dent variables. Note that the rank R can differ from any of the 
ranks of the equations in ( |2.9| ). The rank R must be positive and rational. 

Let V — {ui, W2, . . . ,vq} be the sorted list of all the variables with positive weights, 
including the parameters pi, but excluding ^. The variables are ordered according to 
descending weights: w{vi) is the largest weight, w{vq) is the smallest. The following 
procedure is used to determine the form of the density of rank R : 

(a) Form all monomials of rank R or less by taking combinations of the variables in V. 
Recursively, form sets consisting of ordered pairs. In each pair, the first component 
has a specific combination of different powers of the variables, the second component 
has the weight of the first component. 

Set Bo — {(1; 0)} and proceed as follows: 
For g = 1 through Q do 

For m = through Af — 1 do 

Form Bq^m = M {{Tq^s] Wq_s)}, whcrc M is the number of pairs in Bq-i, 

5=0 

Tq^s = Tq-i^rn Wg^ Wq^s = Wq^i^rn + S w{Vq), (Tq-i^^; Wq^i.m) IS the 

(m + 1)''* ordered pair in Bq-i, and 5q,„ = [ '^,(^^'" 1 
is the maximum allowable power of Vq. 

Set Bq^ [j Bq^ra- 

m=0 

(b) Set Q = Bq. Note that Q has all possible combinations of powers of the variables 
that produce rank R or less. 

(c) Introduce partial derivatives with respect to x. For each pair {Tq s',Wq s) in Q, 

f-PPly '^~i to the term Tq^s, provided £ — R — Wq^s is integer. This introduces just 
enough partial derivatives with respect to x, so that all the pairs retain weight R. 

r)^(T ) 
Gather in a set TC all the terms that result from computing the various — ''^ . 

ox 

(d) Remove those terms from 7i that can be written as a total derivative with respect 

to x, or as a derivative up to terms kept prior in the set. Call the resulting set T, 
which consists of the building blocks of the density p with desired rank R. 

(e) If X has / elements, then their linear combination will produce the polynomial 
density of rank R. Therefore, 

/ 
p = ^CiI{i), 

where T{i) denotes the i*^ element in X, and Ci are numerical coefficients, still to 
be determined. 
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Step 3 Determine the unknown coefficients in the density 

Recall that a conservation law is of the form Dtp + D^J = 0, or Dtp = —D^J, which 
means that Dtp must be the negative of the a:;— derivative of a functional J. 

After computation of Dtp, remove all (ui^t) , j — 0,1,2,... from the expression, 



using the evolution equations in (p.9[). The resulting expression for Dtp must be a total 



X— derivative of some expression. To verify this we use Theorem 2.1 and apply the Euler 
operator. We require that the resulting Euler-Lagrange equations vanish identically by 
the appropriate choice of the coefficients c.^. This leads to a linear system for the q. The 
system must be analyzed and solved for the unknown c.;. In general, the procedure is as 
follows: 

(a) Compute Dtp and replace all [ui^t) ,i = 0, l,...,iV and j = 0,1,2,... using the 



evolution equations in (2.9) 



(b) The resulting expression, called E, must equal Dx{—J) for some functional J. Apply 



to E the Euler operator C from Theorem 2.1. If _E is completely integrable no terms 



will be left, i.e. Cs{E) — 0. If terms remain, set them equal to zero and form the 
linear system for the q, denoted by S. 

(c) Depending on whether or not there are parameters in S, two cases occur: 

(i) If the only unknowns in S are q, solve S for the c^. Substitute the (non-empty) 
solution into the expression of p to obtain its final form. 

(ii) If in addition to the coefficients Ci there are parameters pt in S, then determine 
the conditions on these parameters, so that a density in the given rank exists 
for at least some c^ nonzero. These compatibility conditions assure that the 
system S has non-trivial solutions. Obviously, all Ci equal to zero would give 
the trivial (zero) density. Solving the compatibility conditions may lead to 
different densities of the same rank, corresponding to different choices of the 
parameters. Thus, generating the compatibility conditions enables one to filter 
out all the cases for which there exists a nontrivial density of given rank. Let 
C = {ci, C2, . . . , c/} be the set of all the coefficients that appear in the density. 
In order to determine all possible compatibility conditions, proceed as follows: 

Under the assumption that no pi in S is zero, analyze the system S. First 
determine the coefficients Ci that always must be zero. Exclude these Ci from C 
and set i = i' , where i' is the smallest index of the Cj that remain in C. 

While C^{} do: 

For the building block X{i) with coefficient Ci to appear in p, one needs 
Ci ^ 0. Therefore, set Ci — 1 and eliminate all the other Cj from S. This 
gives compatibility conditions consistent with the presence of the term 
Cil(i) in p. 

If S becomes inconsistent, or compatibility conditions require some of the 
parameters to be zero, 

then: 

Ci must be zero. Hence, set C = C\{ci\, and i — i' , where i' is the 
smallest index of the Cj that remain in C, 
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else: 

Analyze the compatibility conditions and for each resulting branch; 
solve the system S for Cj , and substitute the solution into the expression 
of p to obtain its final form. Then, collect those Cj which are zero for all 
of the branches into a set Z. Since the c^ in Z might not have occurred 
in any of the densities yet, set C — C H Z, and i — i' , where i' is the 
smallest index of the Cj that are still in C. 

(d) Compute J = ^ J E dx, via integration by parts. 

The three examples below illustrate the algorithm. In the first example, we show Steps 1 
and 2 in detail. The second and third example illustrate details of Step 3. The examples 
4.1 and 4.2 illustrate what happens if there are free weights. 



Example 2.2. The wave equation, 

utt - U2x + 'iuu2x + 3uj;^ + auix = 0, (2-10) 

was proposed by Boussinesq to describe surface water waves whose horizontal scale is 
much larger than the depth of the water (Ablowitz and Clarkson, 1991; Hickernell, 1983). 



Conservation laws play a key role in the study of (2.10) for they can be used to prove 
that solutions are bounded for certain sets of initial conditions (Hickernell, 1983), or, 
conversely, to prove that solutions fail to exist after a finite time. 



For computing conservation laws, we rewrite ( ^.10 ) as a system of first-order equations. 



ut + Vx = 0, 
vt + Ux - 3uux - au3x = 0, (2-11) 

where v is an auxiliary dependent variable. It is easy to verify that the terms Ux and au^x 
in the second equation do not allow for uniformity in rank. To circumvent the problem we 
use a trick: we introduce an auxiliary parameter /3 with (unknown) weight, and replace 



(P.llD by 

vt + Pux - 5uux - au3x = 0, (2-12) 



Ut+Vx = 


0, 




3uUx — au^x — 


0, 




Wl,t +U2^x 


= 


0, 


SuiUi^x — OL Ui,3a; 


= 


0. 



or, in our notation 

U2,t + P ui^x-iuiui^x- aui^^x = 0. (2-13) 

Using the procedure in Step 1, we obtain the weights by first computing 

ri,i = lw{^) + lw{ui), ri,2 = l + lw(u2), 

r2,i = lw(f ) + lw(u2), ^2,2 = 1 + 1w(mi) + 1w(/5), 

7-2,3 = 1 + 2w(mi), r2A = 3 + lw(ui), 

then forming the systems 
^1 = {ri,i=ri^2}, M = {?'2,i = ?'2,2 = '"2,3 = ''2,4}, and A = A1UA2, 
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and, finally, solving A for u'(ui), w(u2), w(^) and w{0). This yields 
w{ui) = 2, w(u2) = 3, w(f3) = 2, and w(-—) — 2. 



Hence, the scaling properties of (2.13) are such that 

til ^^ jj ^^ — r- , U2 ^^ TJ q" 1 TTT ^^ "o n" 5 



which expresses that (2.13) is invariant under the scaling symmetry 

{x, t, ui, U2, P) — > (Aa;, X^t, X^^ui, X^^U2, X^"^ (3). 

Let us construct the form of the density with rank i? = 6. Here, V = {u2, ui, /3}, hence, 
111 — U2,V2 — Ml and v^ = (3 and, obviously, Q = 2>. We follow the procedure outlined in 
Step 2: 

(a) For q — \^m — Q: 

&i,o = If] = 2. Thus, with Ti^s — U2'^ , and Wi.s — 3s, where s = 0, 1, 2, we obtain 
ei = Si,o = {(l;0),(u2;3),(u2';6)}. 

For q = 2,m — 0: 

hfi = 1^1 = 3. So, with T2,s = "i", and W2.S = 2s, with s = 0, 1, 2, 3, we obtain 
S2.o-{(l;0),(ui;2),(ui2;4),(ui3;6)}. 

For g = 2, TO = 1: 
we obtain 

^2,1 = {(w2;3),(miM2;5)}, 
since &2,i = [^1 = 1, ^2,5 = M2 ^i^, and W2^s = 3 + 2s, and s = 0, 1. 

For q — 2,m — 2: 

&2,2 = 1^1 - 0. Therefore, ^2,2 = {(m2^;6)}. Hence, 

B2 = {(1; 0), (ui; 2), (ui^; 4), (ui^; 6), (W2; 3), (ui^a; 5), {u2^; 6)}. 

For g = 3: 

we introduce possible powers of (3. So, 

^3.0 = {(l;0),(/3;2),(/32;4),(/33;6)}, B3.4 = {(^2; 3), (/3u2; 5)}, 

B3.1 = {(mi;2),(/3mi;4),(/?V;6)}, ^3,5 = {(^1^2; 5)}, 

S3.2 = {(^i2;4),(/3ui2;6)}, B3.6 = {(^2^6)}. 

^3,3 = {(^1^6)}, 

and 

S3 = {(1; 0), (/?; 2), (/32; 4), (/33; 6), (7/1; 2), (/3wi; 4), (/3^ui; 6), (wi^; 4), 
il3ui^; 6), (ui^; 6), (u2; 3), {(3u2; 5), (^1^2; 5), (1*2'; 6)}. 

(b) Setg = S3. 
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(c) Next, we apply derivatives to the first components of the pairs in Q. 
Computation of £ for each pair of Q gives 

£ = 6,4, 2, 0, 4, 2, 0, 2, 0, 0, 3, 1, 1, and 0. 

Note that in this case aU values for £ are integer. Gathering the terms that come 
from applying the indicated number, £, of partial derivatives with respect to x, gives 

7Y= {0, 13^, UiA^, I3ui^2x, 0^Ui,u\.^,UiUi^2x, I3u\,ul, U2.3x, Pu2,x,UiU2,x,Ui^xU2, uj} . 

(d) Removing from Ti the constant terms, and the terms that can be written as a 
a;— derivative, or as a x— derivative up to terms retained earlier in the set X, yields 

I = {/3^Ui, /?Ui^, Mi^, U2^, Ui,xU2, Ui,x^}. 

(e) Combining these building blocks, the form of the density with rank 6 follows; 

P = Ci /3^Ui + C2 /3mi^ + C3 Ui^ + C4 U2^ + C5 Ui^xU2 + Ce Ui,^^. 

Using Step 3 of the algorithm (illustrated in detail in the next example) we compute the 
density of rank 6 in the original variables: 

p = Pu^ 



Analogously, we computed densities of rank R < 6. The first four densities of ( 2.12 ) are: 



P3 = uv, P4 = Pu 



2 



After substitution of /3 = 1 into the densities above, one gets the densities of (2.11) 
even though initially this system was not uniform in rank. This trick, which involves the 
introduction of one or more extra parameters with weights, can always be attempted if 
any equation in (p.9[) lacks uniformity in rank. 



Example 2.3. Hirota and Satsuma (1981) proposed a coupled system of KdV equations, 

Ut — 6auux + 6vvx — au^^ = 0, 

vt + ^uvx + v3x = 0, (2.14) 

where a is a nonzero parameter. System (p34) describes interactions of two long waves 
with different dispersion relations. It is known to be completely integrable for a = ^. 



In our notation, (2.14) can be rewritten as 

ui,t - 6auiui^x + Gu2U2^x - aui^sx = 0, 

U2,t + iUiU2.x + U2.3X = 0, (2-15) 

which has scaling properties, ui ~ U2 ~ ^tt, -§7 '^ -&^, and a density of rank 4: 

P = Ci Ui^ + C2 UiU2 + C3 U2^- 

Thus far we used Steps 1 and 2 of the algorithm. We now illustrate Step 3, which fixes 
the undetermined coefficients &. 
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(a) Compute Dtp and replace all the mixed derivatives (ui.t) by using (2.15). Then, 

-C2U1 (3uiU2,x + "2,32;) - 2C3M2 {3uiU2,x + U2,3x) ■ 

(b) Apply the Euler operator to get the linear system for the coefficients ci, C2 and C3: 

5 = {(1 + a)c2 = 0, 2 ci + C3 = 0}. 

(c) Obviously, C = {ci,C2,C3} and S has one parameter, a. Thus, we search for com- 
patibility conditions: 

Set ci = 1, which gives 

{Ci = 1, C2 = 0, C3 = -2} 

as one of the solutions without any constraint on the parameter a. Since only 

C2 = 0, one has Z = {02} and C ~ C O Z ~ {C2}, with i' = 2. 

Set C2 = 1. This leads to the compatibility condition a = — 1 and the solution 

{Cl = -- C3, C2 = 1}. 

Since Z — {} and, consequently, C ^ C O Z = {}, the procedure ends. 

Therefore, one gets two densities of rank 4, one without any constraint on a, one with a 
constraint. In summary: 

p = ui^ — 2 W2^, and 

p — — ic3Ui^ + U1U2 + C3U2^, with compatibility condition a — —1. 



A search for densities for (2.14) of rank R < 8 resulted in: 

Rank 2: There is no condition on a. One always has the density p — u. 

Rank 4: At this level, two branches emerge: 

1) Without condition on a, one obtains p — v? — 2u^. 

2) For a = — 1 one has p — uv — ^ {u^ — 2d^), c is free. 

Hence, for a = —1, there is a second independent conserved density: p — uv. 

Rank 6: There is no condition on a. One obtains 

p={l + a)u^ - Zuv'^ - 1^(1 + a.)ux^ + iv^^ ■ 
Rank 8: The system has conserved density 

4 12 22 12 4 2 1 2 " 24 2 " 2 

p = U —U V + -—V - 2uUx + -U2x + -pVUxVx ^WWa; + 71^22; , 

5 5 5 5 5 5 

provided that a—\. 

1 . . 

Therefore, a. = j appears again as one tries to compute the density of rank 8. 



This value of a leads to the only integrable case in the parameterized system (2.14). 
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For a = i, we also computed the density of Rank 10: 



5 20 3 2 ^ 20 4 2 2 ^ 10 2 2 ^ 2 1 2 ^ 40 

p = U ~ —U V + —UV -bU Ux + —V Ux + UU2x - TT"3a; + ^UVM^^fa; 

20 2 2 80 2 2 24 2 4 40 383 

Z-W Wa; - -Z-V Vx ~ —U2xVx - llVU2xV2x + -Z-UV2x - -jVZx ■ 



Example 2.4. Ito (1982) proposed the coupled nonlinear wave equation (Ablowitz and 
Clarkson, 1991): 

Ut + &UUx + 2vVx + U^x = 0, 

vt + 2vux + 2uvx = 0, (2.16) 

which differs from the Hirota-S atsum a system in the interaction and dispersion terms for 

V. 



In the absence of v, system (2.16) reduces to the KdV equation. It is a Hamiltonian 



system with infinitely many conservation laws. The scaling properties of the system ( 2.16 ) 
are 

d^ d d^ 
dx^ ' dt dx^ ' 
and the first five densities are: 
pi — C1U + C2V, p2 = u'^ + v'^, 

P3 = U^+UV^-^Ux^, Pi = U'^+^U^V'^ +\v^ -2uUx^ +\u2x^ -^VUxVx, 

5 324r-2222 2 2 2 2 

P5 = U +—UV +-UV -bU Ux --V Ux +UU2X ^T7"3x - yUTOxWa; -3W Wa; 

2 2 2 

+ ^U2xVx + -jVU2xV2x- 

To illustrate Step 3 of the algorithm in more detail, we consider, 

Ut + &UUx + 2vVx + U3x = 0, 

Vt + a{uxV + uvx) — 0, (2-17) 

which is a parameterized version of (p3£ ) . The form of the density of rank 6 is 



P = Ci U'^ + C2 U^V + C3 UV^ + C4 W"^ + C5 Ux'^ + Cq UxVx + Cy Wj;^. 

We continue with Part (b) of Step 3. After applying the Euler operator, the linear system 
<S for the coefficients ci through C7 is 

S = {(6 — a) C2 = 0, Ci — C3 = 0, 2 C2 — 3 q; C4 = 0, ci + 2 C5 = 0, C3 + 2 C5 = 0, 

cg = 0, C3 + 2 C5 — a C7 = 0, a c-j — Q, 2 C2 + a cg = 0, 2 C2 — 6 cg + a cg = 0}. 

Obviously, C = Ui=i{ci} and S has one parameter, a. Thus, we start the search for 
compatibility conditions. From the sixth and eighth equations in <S we conclude that 
C6 = C7 = 0. Therefore, replace C by C\{cg, cy} ~ lJi;=i{ci}- 

Set ci = 1 and solve 5, which gives 

{Cl = 1, C2 = 0, C3 = 1, C4 == 0, C5 = --, Cg = 0, C7 = 0}, 
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without any constraint on the parameter a. Now, Z — {c2, C4, cg, C7} and C is 

replaced by C fl Z = {c2, C4}, with i' — 2. 

Set C2 — 1. This leads to an inconsistent system. Thus, C ~ C\{c2} = {C4}, with 

z' = 4. 

Set C4 = 1. This implies that a = 0. Hence, C = C\{c4} — {}, and the procedure 

ends. 

In summary: p = u + u v u^ , without any constraint on a. This is the same density 



as for (2.16) 



Computation of the density of rank 8 f or (2.17 ) leads to the condition a — 2 a nd 
density p4 listed above. For a = 2, system ( 2.17 ) reduces to the integrable case ( |2.16| ). 



3. Applications 

For systems with parameters, the algorithm described in Section 2 can be used to find 
the necessary conditions on the parameters so that the systems might have densities of 
fixed rank. If for certain parameters a system has a large number of conserved densities, it 
is candidate to be completely integrable. This is the major application of the algorithm, 
and also of our Mathematica program condens.m. The next examples illustrate the 
computer analysis of equations with parameters. 

3.1. Fifth-Order Korteweg-de Vries Equations 

Consider the family of fifth-order KdV equations, 

ut + au^Ux + (3uxU2x + juusx + u^x = 0, (3.1) 



where a,/?, 7 are nonzero parameters. Special cases of ( |3.lD are well known in the liter- 
ature (Fordy and Gibbons, 1980; Hirota and Ito, 1983; Kupershmidt and Wilson, 1981; 



Satsuma and Kaup, 1977). Indeed, for a = 30, /3 — 20,7 — 10, equation (3.1) reduces to 
the Lax (1968) equation. The SK equation, due to Sawada and Kotera (1974), and also 
Dodd and Gibbon (1977), is obtained for a = 5, /? = 5, 7 = 5. The KK equation, due to 
Kaup (1980) and Kupershmidt, corresponds to a = 20, /5 = 25,7 = 10, a-nd the equation 
due to Ito (1980) arises for a = 2, /3 = 6, 7 = 3. 



The scaling properties of (3.1) are such that 

Q2 Q g5 

dx^ ' dt dx^ 
Using our algorithm, one easily computes the compatibility conditions for the parameters 



a, (3 and 7, so that (3T) admits a polynomial conserved density of fixed rank. 
The results are: 



Rank 2: There are no conditions on the parameters. Indeed, equation (3T) can be 
written as a conservation law with density p ~ u. 

Rank 4: The equation has density p = v? provided that 

/3 = 27. (3.2) 

Only the Lax and Ito equations have this density. 
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Rank 6: The condition 



(2/3^ - 7/37 + 37^) (3.3) 



10 
guarantees the existence of the density of rank 6: 

p = (2/3 - 7)^3 - 15u,2, 7 ^ 2/3. 

The Lax, SK and KK equations have this density. 

Rank 8: There are two branches: 

1) If the condition ( |3.2| ) holds then 

This branch covers the Lax and Ito cases. 

2) If the condition 

a = -^(2/32-7/37-472) (3.4) 

holds, one has 

p = (2/3 + jfu^ - 135(2/3 + j)uux^ + 675^i2x^ 7 ¥^ "2/3. 
This branch covers the SK, KK and Ito equations. 



Rank 10: The conditions 

B = 2t a,nd it = 

10 



/3 = 27 and a= -^7^ (3.5) 



must hold in order to have the density 

p = 7 U — 5O7 U Ux + 1007'U'U2i; U3x 

Only the Lax equation has this density. Naturally, the following que stion arises: What are 



the necessary conditions for the parameters a, (3 and 7 so that (3^) might have infinitely 
many polynomial conservation laws? 

Lax Case : Th e Lax equation admits densities of rank 4 and 6. Combining the conditions 



( P.2D and (^ leads to 

a = ^7' and /3 = 27, (3.6) 

fixing a and (3 in terms of 7 (the latter parameter can always be scaled to 1). The 
conditions in (pj) lead to the Lax equation, which is known to be completely integrable 
and has infinitely many conserved densities. 

SK-KK Cases: The SK and KK equations admit densities of rank 6 and 8. Combining 



(Q and (|J) gives 

1 



a = W and /3 = 7, (3.7) 

5 

a^\^^ and /3 = ^7. (3.8) 

5 2 



The conditions in (3.7) lead to the SK equation, and those in (|3.8|) to the KK equation. 



Both equations are indeed integrable and have an infinite sequence of conserved densities. 
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Ito Case: The Ito equation admits only three densities. Combining (3.2) and (3.4) gives 

2 

a = -7^ and (3 = 27. (3.9) 



These conditions reduce (3.1) to the Ito equation which is not integrable. 
Using our code we were able to retrieve all known integrable cases in the family ( |3.lD , as 
given in e.g. Mikhailov et al. (1987). The Lax equation comes as no surprise since it is the 
fifth-order symmetry of the KdV equation ( p.4| ). Our program confirms that, apart from 
the KdV equation, there are only two non-trivial integrable equations in the family ( |3.lD , 
namely the KK and SK equations. Some of the conserved densities are listed in Tables 
1 and 2. A proof for the gaps (with period three) between conservation laws for the SK 
and KK equations is given in Sanders and Wang (1996). Mathematical rigorous proofs 
for the existence of conservation laws of such evolution equations are given in Sanders 
and Wang (1995b, 1996). 

3.2. Seventh-Order Korteweg-de Vries Equations 

Consider the family of seventh-order KdV equations, 

Ut + au^Uj. + buj' + cuu^U2x + du^u^^ + eu2xU3x + fuxU4x + guu^x + ujx = 0, (3.10) 
where a,b,c,d,e, f and g are nonzero parameters. Special cases of (3.10|) are known in 



t he li terature. For a = 252, b = 63, c = 378, d = 126, e = 63, / = 42,5 = 21, equation 
( ^.10 ) reduces to the SK-Ito equation, due to Sawada and Kotera, and Ito (1980). For 
140, 6 = 70, c = 280, d = 70, e = 70, / = 42, g = 14, equation ( p^ belongs to 



a 



the KdV hierarchy studied by Lax (1968). For a ^ 2016,6 = 630, c = 2268, d = 504, 
e = 252, f = 147, g — 42, the equation belongs to the Kaup-Kupershmidt hierarchy 



(Bilge, 1992; Gerdt, 1993). The scaling properties of ( 3.10 ) are u ~ tt^, tt- ~ tt^- 

ox'' at ox' 

With condens.m, we computed the compatibility conditions for the parameters, so that 

conserved densities of fixed rank exist. The results are: 

Rank 2: p — u \i b — ^[c — 2d). True for both the Lax and SK-Ito equations. 

Rank 4: One has p ^ v? lib = c — id and e — 5{f — 2g); only true for Lax's equation. 

Rank 6: There are three branches. If one of the following three sets of conditions holds: 

a = ^(286/ -42c/ +120/3- 4265 + 63c5-720/2g-fl350/<?2_810<?3), 

5 14c 

d = (2/2-9/<? + 9.g2) and e = 



14^-' '' 2/ -35' 

or 

a = -^(286/ - 42c/- 40/3 -426,g + 63cg + 240/2.g- 450/5^-1-27053), 

d = -l(2/^-9/g + 95-) and , ^ ^ (21c + 20/^ - 9O/5 + 905^) 

42 ^ ■' -ly^ J ' 3 2/ - 35 

or 



a = ^(46/-6c/ + 24d/-665 + 9c5-36d5), 

14c - 14d -I- 10/2 - 45/5 -h 4552 
"" ~ 2733;^ ' 
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then 



Rank 8: 



p = (2/ - 3g)u^ - 2\u^\ f ^ -g. 



p = (49c + 10/^ - 45/5 + '2.Qg^W - 252(2/ - g)uuj + 1764u2a: , 
C7^-^(10/'-45/5 + 20g2), 



provided that 



- — (49c/ + 10/-^ - 196cg - 85/^3 + 2QQfg^ 



^{Uc + 2f-9fg + 4g') 



d^^{7c-2f + 9fg-Ag'), and 



e - 2f-g. 
Combining the conditions for Rank 2 through Rank 8 yields 

a = Y^5^ b = -g'^, c = -g'^, d ^ -g'^, e ^ 3g, f = 2g, 



or 



or 



" " 98^ ' 14^ ' ^ " y^ ' 14^ ' ^ " ^' ■/ = ^5' 



(3.11) 
(3.12) 



4 3-52 9222 « , 7 



(3.13) 
In each case ah of the parameters are fixed in terms of g, which could be scaled to 1. 



The equations (3.11), (3.12), and (3.12) belong to the SK-Ito, Lax and KK hierarchies, 
respectively, and they are all integrable (Bilge, 1992). Some of the conserved densities 
for the SK-Ito and Lax equations are listed in Table 3 for the choices g = 21 and g — 14, 
respectively. 

The conditions in (|3.13| ) are satisfied for a = 2016, b = 630, c = 2268, d = 504, 
e — 252, f = 147,5 = 42, and the first five conserved densities of ( 3.10 ), corresponding 
to the KK-hierarchy, then are: 



Pi 



P2^U^ ~ g^ia;^ 



4 «^ 2 I -'^ 2 

4 4o 



Pi 



i" - —-U-^Ux 



63 



4 2 2 3 

384 "^ 32 ^"^ 1152 ^"^ 



1 



:UU3a; 



7 ""J 4 2 
P5 = U - — U Ma; 



161 



120 640 ^ 2x -r 2ggo 



53 

480 



U U33; 



287 
-320""= 

133 2 1 2 

UOtU^t H UUAt 

5760 240 



97 



192 " 2304 

737 



Uix^, 



UU2x 



1 



17280 



U5x 



Our results merely confirm the computer analysis of (3.10) carried out with REDUCE 
by Bilge (1992) and Gerdt (1993). 



4. More Examples 
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Table 1. Conserved Densities for Sawada-Kotera and Lax 5th-order Equations 



Density 



Sawada-Kotcra equation 



Lax equation 



Pi 



P2 



2" 



P3 



i«3 



i„3 _ 1„ 2 



Pi 



4" 



funaj^ + T_U2x'^ 



TjUUx + onW2i 



PB 



P6 



l„6_«„3„^2_^„^4+6„2„2^ 



i«6 _ |„3„^2 _ 5_^^A _,_ 1„2„2^2 



+2«2i 



i«^ - fu-'ua:^ - fun^;'' + U^U2x'^ 



-I-M8,, 2,,„ 2 , 489 3 _ Mi ,,2„„ 2 , 1 ,, 2,,„ 2 , 10 ,,„„ 3 _ J.,,2,,, 2 

288.., ,,„ 2 I 81,,,,, 2 9 ,,, 2 5 ,,„ ,,„ 2 , 1 „„, , 2 1 „, 2 



^«2a:'U3a:^ + ft««4a:^ " ^«5a:^ 



j«2a:M3a:^ + -^UUix^ - 7^U;,x'^ 



PS 



1„8 _ I,,5,, 2 _ 35,,2„ 4 , 7, ,4 2 



+ |MUa:^«2i:^ + |M^M2a:3 + ^M2a:'' 
+ ^U^UZx'^ - \ux'^UZx^ - ^UU2xU3x^ 
+ j^u'^U4x'^ + ^n2a:M4a:^ - ^UU^x'^ 
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Table 2. Conserved Densities for Kaup-Kupcrshmidt and Ito 5th-order Equations 



Density 



Kaup-Kupcrshniidt equation 



Ito equation 



Pi 



P2 



h' 



P3 



i„3 _ 1„ 2 



Pi 



1„4 _ 9_ 2 I 3_ 2 



i«4 _ |„„,2 + 3, 



pa 



i„6 _ 35 3,, 2 _ .31_,, 4 , 51,,2,,„ 2 



P8 



+ ^«2.3 - ^UUaJ + ^U4J 



P7 



1„7 _ ^„4„^2 _ 369„^^4 + |9„3„2^2 



+ 2619 2,,„ 2 , 2211 3 _ ^2L,,2,,, 2 

+ 4480 "^ "2a: + 2240 ""^a: 1120 ""^a; 



1 71 o 97 o 

-64(j"2-"3a:2 + 36o™4a=2 



-2-«r 2 



i3„6„ 2 _ 427 3,, 4 _ 10431 6 
y ^ - 2 « «a; 32 " "^ 8960 ^^ 



P9 



i„9 _ i3„6„ 2 _ 427 

_21,,5,,„ 2 , 12555 2., 2„„ 2 , 2413 
X 2a; H 448^" "a; «2a; + -324^ 

16461 2„„ 3 , 1641 
+ 1792 "^ "2a; + 356 



_|_ 2413,,3,,„ 3 
+ 224 " "2a; 



IMi „,,„ 4 _ 267 4 2 
256 ""2a ;^i2" '^ 



112" "3^ 

_3699 2 2 _ 4383,,2,,„ ,,„ 2 _ 76635 2 2 
896 ^ 3a; 440 " "2a: "3a: 19712 ^^ 

18891 3 ,141 3,,^ 2 , 8649 2,,^ 2 
19712 "^"3^ + 224" "4^ + 39424"== "4^ 

27639 2 , 2715 3 _ J2;L,,2,,, 2 
+ 19712 ""2^"43: + 39424 "4^ 9856" "5^ 

2943 2 , _9_ 2 9 2 
39424 "2a: "5a -h -^332 ""^^ 39424 "^^ 
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Table 3. Conserved Densities for Sawada-Kotera-Ito and Lax 7th-ordcr Equations 



Density 



Sawada-Kotera-Ito equation 



Lax equation 



Pi 



P2 



i«2 



P3 



1„3 _ 1,, 2 



3 6 



Pi 



4" 4 



««i: + T9"2i 



P5 



iw"' - M^Ml^ + j-UU2x^ ~ ■fo'^3^^ 



P6 



PT 



i«^ - ■iu*Ux^ - f M«,4 + i|„3„2,2 



i«^ - lu^uj 



_l_24„ 2,,„ 2 , 163, ,„o 3 _ -29.„2,,,, 2 , 1„ 2„„ 2 _, 10,,„„ 3 _ J.,,2,,,, 

35 31^ 2a: — T05 ^a? ^2 ^a: ^ 21 2a; — 14 3a: 



32 9 1 9 1 9 5 9 1 

-^«2a;M3a; + 35«"4i; - gJgMSa; - -^U2xU3x + -^UUAx 



r«5a 



P8 



1 ,,8 _ I„5,, 2 _ 35„2,, 4 , I„4,,„ 2 



-^UUx'^U2x'^ + |M^W2a;3 + ^"2a;'' 
-^U^USx'^ - jUx'^USx'^ - ^UU2xU3x^ 
--^U^Uix^ + j^U2xU4x'^ - -i^UUSx'^ 



yUe,x 
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4.1. Nonlinear Schrodinger Equation 

The nonlinear Schrodinger (NLS) equation (Ablowitz and Clarkson, 1991), 

iqt-q2x + 2\qfq^0, (4.1) 

arises as an asymptotic Hmit of a slowly varying dispersive wave envelope in a nonlinear 
medium, and as such has significant applications in nonlinear optics, water waves, and 



plasma physics. Equation (4.1) is known to be completely integrable, and together with 



the ubiquitous KdV equation (UM, is one of the most studied soliton equations. 



There are two different ways to compute conserved densities for (4.1). 
Method 1 One can replace (^J) by 

Ut—V2x+'^v{u +v) ~ 0, 

vt + U2x - Mu^ + y'^) = 0, (4.2) 

where q — u + iv. The scaling properties are such that 

d d d^ 
dx ' dt dx^ 
We computed seven conserved densities with our program, of which the first six are: 



/35 = U^ + Zu^V^ + Zu^V + V^ + bU^Ux^ + iv^Ux^ + 'iv^Vx^ + hv^Vx^ + 4 UVUxVx 



+ ^U2x + 2*^20. , 




























3 


+ 


1 5 

5 


+ 


1 

-;:UUx 
3 


9 

Vx 


+ 


1 

3" 


'^U2xVx 


+ 


1 2 

o^ U2xVx 


+ 


1 

3 


+ ^^U,xV2.. 





























Pq = U VUx + TtU V Ux + —VUx + -^UUx Vx + —UU2xVx + —VU2xVx + ■7:VUxVx 



Method 2 One could consider q and q* as independent variables and add the complex 



conjugate of (4.1) to the NLS equation. Absorbing i in the scale of t then yields: 

qt ~ q2x + '2q^q* = 0, 

q:+q*2x-2q*^q = 0. (4.3) 

According to the procedure described at the end of Step 1 in Section 2.3, our program 
computes the constraints 

wiq) = 2-wiq*), w{q*)^w{q*), ^(^) = 2, 

and sets the left hand sides of the first two equal to one. Then, it solves the equations 

l = 2-wiq*), l^wiq*), 

piece by piece. Both lead to the solution w{q*) = 1. Hence, the program continues with 
w{q) ^w{q*) = 1. 
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For rank 2 through rank 6 our program produces the conserved densities: 
Pi = m*, P2 = 9*fe, 

P3 = q'^q*'^ + QxqI, pa = qq*'^qx + \q2xql, 

3*3|-'^*2 2|^ * *i^2*2|^ * 

P5 = q q +2*? 9x + 4<?<? g^g^r + ^q qx + 2q2xq2x- 

Obviously, these two sets of conservation laws are connected (but not piece by piece) via 
a simple change of variables: -u = 5(9 + 9*), t' = -^{q ^ q*)- 

The second method has the advantage that the conserved densities are expressed in 
the original variable q and it conjugate q* . On the other hand, the conserved densities 
from Method 2 may not be independent. 

4.2. Non-dispersive Long Wave System 

The non-dispersive long wave equations (Kupershmidt, 1985a) 

Ut + VUx + UVx = 0, 

Vt+Ua:+ VV^ = 0, (4.4) 

is another example of an integrable system. 

We use this example to illustrate how our code determines a free weight. Indeed, for 



(4.4), w{u) = 2w{v), w{v) — w{v), w{^) = 1-1- w(v), with w{v) as the only free weight. 
As described in Step 1 of the algorithm, the program sorts the right hand sides of these 
constraints, sets their left hand sides equal to one, and proceeds with solving 

1 — w{v), 1 = 2w{v), 1 = w{v) + 1, 

one by one. That leads to the choices w{v) = 1, w{v) = i, and wlv) = 0. Since the latter 
is zero it is incremented by one to get w{v) — 1. 

The first choice, w{v) = 1, gives w{u) = w{-^) = 2. The second choice, w{v) — \, gives 
«;('«) = l,w(^) = |. Another valid choice (not considered by our program) would be 
w{v) — J, wiu) — 27 and w{-^) — J. Obviously, there are infinitely many fractional 
choices for w{v). Recall that fractional weights and ranks are indeed allowed. 

The program continues automatically with the smallest integer choice, w{v) = 1. 
Hence, w{u) = 2, or, symbolically, 

ax^ ox 

For rank one through eight we obtained the following densities: 

P2 = M, 

P8 = u"^ + Gm^w^ -I- iu^v^ \- \uv^ . 

This set of densities remains the same for any valid choice oiw(v). Indeed, we could have 
computed the conserved density p — uv with different choices of the free weight w(y). 
This p has rank _R — w{u) + w{v) = 3w{v). If we choose w{v) = 2 then p = uv has 
rank i? = |. On the other hand, for the choice w{v) — ^, this density has rank R — j. 



Pi 


z:^ 


V, 


PS. 


= 


uv, 


P5 


= 


u^v -t- ^uv^ , 


P7 


= 


u^v + u^v'^ + 11)""^ 
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In conclusion, choosing w{v) differently does not affect the form of the desired densities, 
provided the rank of p is also adjusted appropriately. Due to the superposition principle, 
the same argument can be made if p is the sum of many terms, involving x— derivatives 
of the dependent variables. 

4.3. Three-Component Korteweg-de Vries Equation 

Consider the 3-component extension of the KdV equation (Kupershmidt, 1985b), 

Ut — 6uUx + 2vVx + 2wWx — U3x = 0, 

vt — 2vUx — 2uVx = 0, (4-5) 

Wt — 2wUx — 2uWx — 0, 

which can be written as a bi-Hamiltonian system with infinitely many conservation laws. 
The scaling properties indicate that 

U ^^^ V ^^^ IV ^^^ ^"^ 

dx"^ ' dt dx^ ' 
and the first four densities for (|4.5|) are: 



Pl = ClU + C2V + C-^W , P2 — U ^V — W , P3 — U — UV ~ UW Ux^ ^ 

4 ^22,4 ^22,22,4 o 2,2, , 

Pa = U U %] -\-—V U W +—V W +— W —lUUrr +—Ur,^ + —VUT.VT + —WUxWT.. 

^ 5 55 5 5 X ^ 2x ^ X X ^ X X 

Obviously, from pi we can see that u, v and w are independent conserved densities. 

5. Using the Program condens.m 

We now describe the features, scope and limitations of our program condens.m, which 
is written in Mathematica syntax (Wolfram, 1996). The program condens.m has its 
own menu interface with 30 samples of data files. Users are assumed to have access to 
Mathematica (version 2.0 or higher). The code condens.m and the data files should be 
put in one directory. 

5.1. The Menu Interface 

After Mathematica comes up with 'In[l]:=', type 

In[l]:= <<condens.m 

to read in the code condens.m and start the program. Via its menu interface, the 
program will automatically prompt you for answers. 

Example 5.1. Let us compute the density of rank 4 (provided it exists) for the Drinfel'd- 
Sokolov system (Ablowitz and Clarkson, 1991): 

Ut + 3vvx = 0, 

Vt + 2v3x + 2uvx + UxV ^ 0. (5.1) 

Since this example is in the menu, start the program and pick entry 25 from the menu. 
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*** MENU INTERFACE *** (page: 3) 

21) MVDNLS System (d_mvdiils .m) 

22) Kaup System-parameterized (d_pkaup.m) 

23) Kaup System (d_kaup.m) 

24) Kaup-Broer System (d_broer.m) 

25) Drinf el'd-Sokolov System (d_soko.m) 

26) Dispersive Long Wave System (d_disper.m) 

27) Non-dispersive Long Wave System (d_nodisp.m) 

28) 3-Component KdV System (d_3ckdv.m) 

29) 2-Component Nonlinear Schrodinger Equation (d_2cnls.m) 

30) Boussinesq System (d_bous.m) 
nn) Next Page 

tt) Your System 

qq) Exit the Program 

ENTER YOUR CHOICE: 25 

Enter the rank of rho : 4 

Enter the name of the output file: drisokr4.o 

WELCOME TO THE MATHEMATICA PROGRAM 
by UNAL GOKTAS and WILLY HEREMAN 
FOR THE COMPUTATION OF CONSERVED DENSITIES. 
Version 3.0 released on February 24, 1997 
Copyright 1996 



2 

The normalized density rho[x,t] is : u 

2 

2 2 

The corresponding flux j[x,t]is: 2u u -2u +4u u 

1 2 2,x 2 2, XX 

Result of explicit verification (rho_t + J_x) = 
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At the end of computation, the norniahzed density and flux are available. In the absence 
of parameters, both are normalized according to the coefficient of the first term in the 
density. If there are parameters, the common denominators that have been multiplied 
through are shown. To see the density and flux in standard Mathematica notation, type 
rho[x,t] or j[x,t]. However, type pdef orm[rho [x,t] ] or pdef orm[j [x,t] ] to see 
them in subscript notation. 

Note that the form of the densities p is not unique. Densities can always be integrated 
by parts to obtain equivalent forms, modulo total derivatives. In Mathematica version 
2.2, equivalent forms can be obtained via the command Integrate [rho[x,t] ,x]. 

5.2. Preparing Data Files 

To test systems that are not in the menu, prepare a data file in the format of our data 
files. Of course, the name for a new data file should not coincide with any name already 
listed in the menu, unless you intended to modify those data files. 



Example 5.2. For the parameterized Hirota-Satsuma system (2.14) the data file reads: 

(* start of data file d_phrsat.m *) 

debug = False; 

(* Hirota-Satsuma System with one paraimeter *) 

eq[l] [x,t]=D[u[l] [x,t] ,t]-aa*D[u[l] [x,t] ,{x,3}]- 

6*aa*u[l] [x,t]*D[u[l] [x,t] ,x]+6*u[2] [x,t]*D[u[2] [x,t] ,x] ; 

eq[2] [x,t]=D[u[2] [x,t],t]+D[u[2] [x,t] ,{x,3}] +3*u[l] [x,t]*D[u[2] [x,t] ,x] ; 

noeqs = 2; 

name = "Hirota-Satsuma System (parameterized)"; 

parameters = {aa}; 

weightpars = {}; 

(* user can supply the rank of rho and a neune for the output file *) 

(* rhorank =6; *) 

(* myfile = "phrsatr6.o; *) 

(* user can give the weights of u[l] , u[2] and partial t *) 
(* weightu[l]=2; weightu[2] =2; weight [t] =3; *) 

(* user can give the form of rho. Here, for density of rank 6: *) 
(* formrho[x,t]={c[l]*u[l] [x,t] '3+c [2] *u[l] [x,t]*u[2] [x,t] ~2+ 

c[3]*D[u[l] [x,t],x]-2+c[4]*D[u[2] [x,t],x]-2}; *) 

formrho[x,t] = {}; 

(* end of data file d_phrsat.m *) 
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Explanation of the lines in the data file: 

debug = False; 

No trace of internal computations, otherwise, set it True. 

eq[k] [x,t] = . . . ; 

Give the k^^ equation of the system in Mathematica notation. Note that there is no == 
at the end. 

noeqs = 2; 

Specify the number of equations in the system. 

name = "Hirota-Satsuma System (parameterized)"; 
Pick a short name for the system. The quotes are essential. 

parameters = {aa}; 

Give the list of the parameters in the system. If there are none, set parameters ~ { }■ 

weightpars = {}; 

Give the list of the parameters that are assumed to have weights. Note that weighted 
parameters are not listed in parameters, which is the list of parameters without weight. 

(* rhorank =6; *) 

Optional; give the desired rank of the density. Useful if you want to work with the 
program less interactively (in batch mode). 

(* myfile = "hirsatr6.o; *) 

Optional; give a name of the output file. Useful for less interactive use of the program. 

(* weightu[l]=2; weightu[2] =2; weight [t] =3; *) 

Optional; specify a choice for some or all of the weights. The program then skips the 

computation of the weights, but still checks for consistency. Particularly useful if there 

are several free weights and non-interactive use is preferred. 

f ormrho [x,t] = {}; 

The program will compute the form of p of the given rank. 

formrho[x,t]={c[l]*u[l] [x,t] "3+c [2] *u[l] [x,t]*u[2] [x,t]'2+ 
c [3] *D [u [1] [x , t] , x] -2+c [4] *D [u [2] [x , t] , x] "2} ; 

Alternatively, one could give a form of p (here for rank 6). The density must be given 
in expanded form and with coefficients c[i]. The braces are essential. If p is given, the 
program skips both the computation of the weights and the form of the density. Instead, 
the code uses what is given and computes the coefficients c[i]. This option allows one, for 
example, to test densities from the literature. 

Anything within (* and *) are comments, and ignored by Mathematica. 

Once the data file is ready, run it via the choice "tt) Your System" in the menu. 
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5.3. Scope and Limitations 

Our program can handle PDEs that can be written as systems of evolution equations. 
The evolution equations must be polynomials in the dependent variables (no integral 
terms). Only two independent variables {x and t) are allowed. No terms in the evolution 
equations should explicitly depend on x and t. 

Theoretically, there is no limit on the number of evolution equations. In practice, for 
large systems, the computations may take a long time or require a lot of memory. The 
computational speed depends primarily on the amount of memory. 

The program only computes polynomial conserved densities in the dependent variables 
and their derivatives, without explicit dependencies on x and/or t. 

By design, the program constructs only densities that are uniform in rank. The uniform 
rank assumption for the monomials in p allows one to compute independent conserved 
densities piece by piece, without having to split linear combinations of conserved densities. 
Due to the superposition principle, a linear combination of conserved densities of unequal 
rank is still a conserved density. This situation arises frequently when parameters with 
weight are introduced in the PDEs. 

The input systems may have one or more parameters, which are assumed to be nonzero. 
If a system has parameters, the program will attempt to compute the compatibility 
conditions for these parameters such that densities (of a given rank) exist. 

The assumption that all parameters in the given evolution equation must be nonzero 
is necessary. As a result of setting parameters to zero in a given system of evolution 
equations, the weights and therefore the rank of p might change. 

In general, the compatibility conditions for the parameters could be highly nonlinear, 
and there is no general algorithm to solve them. The program automatically generates 
the compatibility conditions, and solves them for parameters that occur linearly (see 
Section 3.2). Grobner basis techniques could be used to reduce complicated nonlinear 
systems into equivalent, yet simpler, non-linear systems. For PDEs with parameters and 
when the system for the coefficients Ci is complicated, the program saves that system 
and its coefficient matrix, etc., in the file worklog.ni. Independent from the program, the 
worklog files can later be analyzed with Mathematica functions. 

The assumption that the evolution equations are uniform in rank is not very restrictive. 
If the uniform rank condition is violated, the user can introduce one or more parameters 
with weights. This also allows for some fiexibility in the form of the densities. Although 
built up with terms that are uniform in rank, the densities do not have to be uniform in 
rank with respect to the dependent variables alone. This is illustrated in Example 2.2. 

Conversely, when the uniform rank condition is fulfilled, the introduction of extra 
parameters (with weights) in the given PDE leads to a linear combination of conservation 
laws, not to new ones. 

In cases where it is not clear whether or not parameters with weight should be intro- 
duced, one should start with parameters without weight. If this causes incompatibilities 
in the assignment of weights (due to non- uniformity) , the program may provide a sug- 
gestion. Quite often, it recommends that one or more parameters be moved from the list 
of parameters into the list weightpars of weighted parameters. 

For systems with two or more free weights, the user will be prompted to enter values for 
the free weights. If only one weight is free, the program will automatically compute some 
choices for the free weight, and eventually continue with the lowest integer or fractional 
value (see Examples 4.1 and 4.2). The program selects this value for the free weight; it 
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is just one choice out of possibly infinitely many, ff the algorithm fails to determine a 
suitable value, the user will be prompted to enter a value for the free weight. 

Negative weights are not allowed, except for w{-^). If w{-^) < 0, the program will 
give a warning, but continue with the computations. Zero weights are allowed, but at 
least one of the dependent variables must have positive weight. The code checks these 
conditions, and if they are violated the computations are aborted. Note that fractional 
weights and densities of fractional rank are permitted. 

Our program is a tool in the search of the first half-dozen conservation laws. An 
existence proof (showing that there are indeed infinitely many conservation laws) must be 
done analytically, e.g. by explicitly constructing the recursion operator (Kodama, 1985; 
Sanders and Wang, 1995b, 1996) that connects conserved densities, or by computing 
high-order symmetries with Lie symmetry software (Hereman, 1996). If our program 
succeeds finding a large set of independent conservation laws, there is a good chance 
that the system of PDEs has infinitely many conserved densities and that the recursion 
operator could be constructed explicitly. If the number of conservation laws is 3 or less, 
most likely the PDEs are not integrable, at least not in that coordinate representation. 



6. Other Software Packages 

This section gives a review of other symbolic software for the computation of conserved 
densities, together with recent developments in this area. 

6.1. SYMCD 

The program SYMCD, written by Ito (1994), is an improved version of CONSD by 
Ito and Kako (1985). Both programs are in REDUCE. 

Similar to our program, SYMCD uses scaling properties to compute polynomial con- 
served densities for systems of any number of evolution equations with uniform rank. 
CONSD had a limit on the number of evolution equations. This limitation was removed 
in SYMCD. Evolution equations must be polynomial in the dependent variables and 
their derivatives, and variables with negative weight are not allowed. With REDUCE 3.5 
on an IBM Rise 6000 work station, we tested the version of SYMCD released in 1996. 
Our test cases included equations (p!^, ( ^.12[ ), ( |2.14| ), (|j]) and (^^. 



For systems with or without parameters, SYMCD gives the same conserved densities 
as our program (up to terms that can be converted via integration by parts). 

However, SYMCD does not properly handle systems with parameters. It stops after 
generating the necessary conditions on the parameters, which must be analyzed sepa- 
rately. Such analyses revealed that the conditions do not always lead to a density of fixed 
rank. Indeed, in solving for the undetermined coefficients, SYMCD considers all possible 
branches in the solution, irrespective of whether or not these branches lead to a conserved 
density, as confirmed by Ito (1996). Another major difference is that parameters with 
weights are not allowed in SYMCD, which restricts the scope of SYMCD to systems 
with uniform rank. In conclusion, our code is more sophisticated than SYMCD in han- 
dling systems with parameters and systems of non-uniform rank. For more information 
contact Ito at ito@puramis.amath.hiroshima-u.ac.jp. 
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6.2. DELiA 



The PC package DELiA for Differential Equations with Lie Approach^ developed in 
the period 1989-1991, is an outgrowth of the SCoLar project by Bocharov and Bronstein. 
DELiA, written in Turbo PASCAL by Bocharov (1991) and coworkers, is a commercial 
computer algebra system for investigating differential equations using Lie's approach. 

The program deals with symmetries, conservation laws, integrability and equivalence 
problems. It has a special routine for systems of evolution equations, which we used for 
computing conserved densities. We tested DELiA 1.5.1 on a P C wit h Intel 48 6 pro cessor 
and 16 MB of RAM. Our test cases included equations dJ), ( ^.12[ ), ( |2.14| ), (|2.14| ) with 
a = 1/ 2, (O ) with arbitrary a, /? = 20 and 7 = 10, (Q and (|4.5|). 

For (2^) and (|3j) with a = 1/2, DELiA returned the same densities as the ones 
listed in this paper, up to ter ms th at differ via i nteg ration by parts. With DELiA, one 
cannot compute densities for ( 2.12 ), (Lj), and (L5). These systems are out of its class: 
the program requires second or higher-order spatial derivative terms in all equations. 
For systems with parameters, DELiA does not automatically compute the densities 
corresponding to the (necessary) conditions on the parameters. One has to use DELiA's 
integrability test first, which determines conditions based on the existence of formal 
symmetries. Since these integrability conditions are neither necessary nor sufficient for 
the existence of conserved densities, one must further analyze the conditions manually. 
Once the parameters are fixed, one can compute the densities. For further information 
we refer to Bocharov at alexei@wri.com. 



6.3. FS 



The REDUCE program FS for "formal symmetries" was written by Gerdt and Zharkov 
(1990). The code FS can be applied to polynomial nonlinear PDEs of evolution type, 
which are linear with respect to the highest-order spatial derivatives and with non- 
degenerated, diagonal coefficient matrix for the highest derivatives. The algorithm in FS 
requires that the evolution equation are of order two or higher in the spatial variable. 
However, the formal symmetry approach does not require that the evolution equations 
are uniform in rank. 

We teste d FS with REDUCE 3.5 on an IBM Rise 6000 work station for equations 
(2.4), ( 2.12|) , (3T), (4.4) and (L5). We were unable to compute the densiti es f or systems 
( 2.12[ ), (4.4), and (15), since FS is not apphcable to such systems. For (^.4|), FS gave 
the same densities as our program, up to terms that differ via integration by parts. 

Like DELiA, applied to equations with parameters, FS computes the conditions on the 
parameters using the symmetry approach. For ( |3.l| ), FS and condens.m gave equivalent 
sets of compatibility conditions on the parameters. More information can be obtained 
from Gerdt at gerdt@jinr.dubna.su. 



6.4. Software under Development 



In addition to the available software, several research groups are currently developing 
software for conserved densities. 

Sanders and Wang at the Free University of Amsterdam are developing a software 
package in Maple and FORM to compute both conserved densities and recursion oper- 
ators (Sanders and Roelofs, 1994; Sanders and Wang, 1995a). Their approach is more 
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abstract and relics on the implicit function theorem. In contrast to our algorithm, their 
code makes no assumptions about the form of the conserved density. Further information 
can be obtained from Sanders at jansa@cs.vu.nl. 

Wolf (1995) at the University of London is designing a package in REDUCE for the 
computation of conserved densities, which implements two methods. If the density is of 
fixed order, then the program constructs it directly. This method is currently limited 
to two independent variables. In the second method the characteristic function of the 
conservation law is constructed first, and the density and flux are recovered via integra- 
tion by parts. There is no limitation on the number of independent variables for this 
method. Both approaches use Wolf's program CRACK for solving overdetermined sys- 
tems of PDEs. See Hereman (1996) for a review of CRACK and for additional references 
to Wolf's work. Wolf's algorithm is particularly efficient for showing the non-existence 
of conservation laws of high order. In contrast to our program, it also allows one to 
compute non-polynomial conserved densities. For further information contact Wolf at 
T.Wolf@maths.qmw.ac.uk. 

Ahner, Tschantz, and Cook at Vanderbilt University are working on a similar project in 
Mathematica (Ahner, 1995). We have no details about their algorithm or implementation. 
Contact Ahner at ahnerjf@ctrvax.vanderbilt.edu for further information. 

Hickman (1996) at the University of Canterbury, Christchurch, New Zealand, has im- 
plemented a slight variation of our algorithm in Maple. Instead of computing the differ- 
ential monomials in the density by repeated differentiation, Hickman uses a tree struc- 
ture combining the appropriately weighted building blocks. For further information, send 
email to M.Hickman@math.canterbury.ac.nz. 

7. Conclusions 

A prototype of condens.m was used to compute conserved densities of a generalized 
KdV equation due to Schamel (Verheest and Hereman, 1995), and a modified vector 
derivative NLS equation (Willox et ai, 1995). Based on the results obtained with the 
software, integrability questions for these equations could be answered. 

We offer the scientific community a Mathematica package to carry out the tedious 
calculations of conserved densities for systems of nonlinear evolution equations. Details 
about the algorithm and its implementation can be found in G6kta§ (1996a), and G6kta§ 
and Hereman (1997). The code condens.m, together with several data and output files, 
is available via anonymous FTP (G6kta§, 1996b). 

Extensions of the algorithm presented in this paper towards PDEs with more than one 
spatial variable, dynamical systems, and systems of differential-difference equations are 
under development. Further information about extensions can be obtained via email to 
ugoktas@mines.edu or whereman@mines.edu. 
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